* This file carries out the analysis in Table 2 (demographic tests) and outputs it

use "$saveddata/data_complete_withcosts.dta", clear 

set more off




		preserve	
		* Collapse to correct level
		collapse (count) freq = aeattenddisp (mean) d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat, by(bin)

		* Now do analysis...
		egen sumfreq = sum(freq)
		gen pfreq = freq/sumfreq
		drop sumfreq

		lab var pfreq "Proportion of patients"
		lab var bin "Exit time (10 min intervals)"
		compress

		** ANALYSIS HERE
		
		* bin polynomial variables
	forv x = 1/15 {
		gen bin`x' = bin^`x'
		}

	* parameters
	
	local porder = 10
	local start = 180
	local end = 240

	* waiting times - find end point
	local sumdiff = 999
	while abs(`sumdiff')>0.01 {
		
		cap drop pfreq_cf diff
		
		* display then update end point
		di "End point: `end'"
		local end = `end' + 10
		
		* compute counterfactual
		qui reg pfreq bin1-bin`porder' if bin<`start' | bin>`end'
		predict pfreq_cf, xb

		* check difference in mass
		gen diff = pfreq - pfreq_cf 
		qui su diff if bin>`start' | bin<`end'
		local sumdiff = r(sum)
			
		
		}
	drop diff

	* waiting times - compute effects
	qui reg freq bin1-bin`porder' if bin<`start' | bin>`end'
	predict freq_cf, xb
	gen wait = freq*bin 

	gen wait_cf = freq_cf*bin 
	gen wait_diff = wait_cf - wait

	su wait_diff
	local saving = r(sum)
	su freq
	local patients = r(sum)/2
	su wait_cf
	local totalwait = r(sum)/2
	su wait_diff if bin>=`start' & bin<=`end'
	local saving_int = r(sum)
	su wait_cf if bin>=`start' & bin<=`end'
	local totalwait_int = r(sum)/2
	su freq if bin>=`start' & bin<=`end'
	local patients_int = r(sum)/2
	local savingpp = round(`saving'/`patients',1)
	local savingpc = round(`saving'/`totalwait',0.01)
	local savingpp_int = round(`saving_int'/`patients_int',1)
	local savingpc_int = round(`saving_int'/`totalwait_int',0.01)
	di "Total waiting time saving (mins per patient): " `savingpp'
	di "Total waiting time saving (% of c.f. total): " `savingpc'
	di "Total waiting time saving in distorted window (mins per patient): " `savingpp_int'
	di "Total waiting time saving in disorted window (% of c.f. total): " `savingpc_int'

	* weights for pre-threshold test
	egen sumfreq_pre_obs = sum(freq) if bin>`start' & bin<=240
	gen w_pre_obs = freq/sumfreq_pre_obs

	egen sumfreq_pre_cf = sum(freq_cf) if bin>`start' & bin<=240
	gen w_pre_cf = freq_cf/sumfreq_pre_cf

	gen freq_diff = freq_cf - freq 
	egen sumfreq_post_diff = sum(freq_diff) if bin>240 & bin<=`end'
	gen w_post_diff = freq_diff/sumfreq_post_diff

	gen nonmover_pre = sumfreq_pre_cf/sumfreq_pre_obs
	
	* weights for DP test
	egen sumfreq_dp_obs = sum(freq) if bin>`start'& bin<=`end'
	gen w_dp_obs = freq/sumfreq_dp_obs

	egen sumfreq_dp_cf = sum(freq_cf) if bin>`start'& bin<=`end'
	gen w_dp_cf = freq_cf/sumfreq_dp_cf

	* other outcomes - compute counterfactuals and selection-only counterfactuals
	foreach var of varlist d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat {
	


		* counterfactual outcomes
		
		local porder = 10
		
		qui reg `var' bin1-bin`porder' if bin<`start' | bin>=260
		predict `var'_cf, xb
		
		* pre-threshold test

			* non-movers: pre-target counterfactual average outcome 
			qui egen pre_`var'_pre_cf = sum(`var'_cf*w_pre_cf)
				
			* movers: post-target counterfactual average outcome (weighted by mover proportions)
			qui egen pre_`var'_post_cf = sum(`var'_cf*w_post_diff)

			* composition-adjusted counterfactual
			qui gen pre_`var'_cac = nonmover_pre*pre_`var'_pre_cf + (1-nonmover_pre)*pre_`var'_post_cf
			
			* observed outcome 
			qui egen pre_`var'_obs = sum(`var'*w_pre_obs)
			
		* DP test
	
			* counterfactual
			qui egen dp_`var'_cf = sum(`var'_cf*w_dp_cf)
				
			* observed outcome 
			qui egen dp_`var'_obs = sum(`var'*w_dp_obs)

		}


	format %10.0fc freq* sumfreq*
		
	* build results table
	gen helper = 1
	collapse (mean) nonmover_pre pre_* dp_*, by(helper)

	foreach i in  d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat {
		
		gen diff_`i' = pre_`i'_obs - pre_`i'_cac
		gen diff_pc_`i' = diff_`i' / pre_`i'_cac
		
		gen dp_diff_`i' = dp_`i'_obs - dp_`i'_cf
		gen dp_diff_pc_`i' = dp_diff_`i' / dp_`i'_cf
		}
		
	foreach x in d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat {

		su diff_`x'
		matrix observe_`x' = r(mean)
		su diff_pc_`x'
		matrix observe_pc_`x' = r(mean)
		
		su dp_diff_`x'
		matrix observe_dp_`x' = r(mean)
		su dp_diff_pc_`x'
		matrix observe_pc_dp_`x' = r(mean)
		}
		

restore

set more off

capture program drop myboot

program define myboot, rclass

* data for aggregate analysis
	preserve

	* draw bootstrap sample
	gen trust_code = substr(site_code,1,3)
	bsample, cluster(trust_code)
	
	cap drop diff_death
	
set more off
		* Collapse to correct level
		collapse (count) freq = aeattenddisp (mean) d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat, by(bin)

		* Now do analysis...
		egen sumfreq = sum(freq)
		gen pfreq = freq/sumfreq
		drop sumfreq

		lab var pfreq "Proportion of patients"
		lab var bin "Exit time (10 min intervals)"
		compress

		** ANALYSIS HERE
		
		* bin polynomial variables
	forv x = 1/15 {
		gen bin`x' = bin^`x'
		}

	* parameters
	local porder = 10
	local start = 180
	local end = 240

	* waiting times - find end point
	local sumdiff = 999
	local t = 0
	while abs(`sumdiff')>0.01 & `t'<100 {
		
		cap drop pfreq_cf diff
		
		* display then update end point
		di "End point: `end'"
		local end = `end' + 10
		
		* compute counterfactual
		qui reg pfreq bin1-bin`porder' if bin<`start' | bin>`end'
		predict pfreq_cf, xb

		* check difference in mass
		gen diff = pfreq - pfreq_cf 
		qui su diff if bin>`start' | bin<`end'
		local sumdiff = r(sum)
			
		local t=`t'+1
		}
	drop diff

	* waiting times - compute effects
	qui reg freq bin1-bin`porder' if bin<`start' | bin>`end'
	predict freq_cf, xb
	gen wait = freq*bin 

	gen wait_cf = freq_cf*bin 
	gen wait_diff = wait_cf - wait

	su wait_diff
	local saving = r(sum)
	su freq
	local patients = r(sum)/2
	su wait_cf
	local totalwait = r(sum)/2
	su wait_diff if bin>=`start' & bin<=`end'
	local saving_int = r(sum)
	su wait_cf if bin>=`start' & bin<=`end'
	local totalwait_int = r(sum)/2
	su freq if bin>=`start' & bin<=`end'
	local patients_int = r(sum)/2
	local savingpp = round(`saving'/`patients',1)
	local savingpc = round(`saving'/`totalwait',0.01)
	local savingpp_int = round(`saving_int'/`patients_int',1)
	local savingpc_int = round(`saving_int'/`totalwait_int',0.01)
	di "Total waiting time saving (mins per patient): " `savingpp'
	di "Total waiting time saving (% of c.f. total): " `savingpc'
	di "Total waiting time saving in distorted window (mins per patient): " `savingpp_int'
	di "Total waiting time saving in disorted window (% of c.f. total): " `savingpc_int'

	* weights for pre-threshold test
	egen sumfreq_pre_obs = sum(freq) if bin>`start' & bin<=240
	gen w_pre_obs = freq/sumfreq_pre_obs

	egen sumfreq_pre_cf = sum(freq_cf) if bin>`start' & bin<=240
	gen w_pre_cf = freq_cf/sumfreq_pre_cf

	gen freq_diff = freq_cf - freq 
	egen sumfreq_post_diff = sum(freq_diff) if bin>240 & bin<=`end'
	gen w_post_diff = freq_diff/sumfreq_post_diff

	gen nonmover_pre = sumfreq_pre_cf/sumfreq_pre_obs
	
	* weights for DP test
	egen sumfreq_dp_obs = sum(freq) if bin>`start'& bin<=`end'
	gen w_dp_obs = freq/sumfreq_dp_obs

	egen sumfreq_dp_cf = sum(freq_cf) if bin>`start'& bin<=`end'
	gen w_dp_cf = freq_cf/sumfreq_dp_cf

	* other outcomes - compute counterfactuals and selection-only counterfactuals
	foreach var of varlist d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat{
	
		* counterfactual outcomes
	
		local porder = 10
		
		qui reg `var' bin1-bin`porder' if bin<`start' | bin>=260
		predict `var'_cf, xb
		
		* pre-threshold test

			* non-movers: pre-target counterfactual average outcome 
			qui egen pre_`var'_pre_cf = sum(`var'_cf*w_pre_cf)
				
			* movers: post-target counterfactual average outcome (weighted by mover proportions)
			qui egen pre_`var'_post_cf = sum(`var'_cf*w_post_diff)

			* composition-adjusted counterfactual
			qui gen pre_`var'_cac = nonmover_pre*pre_`var'_pre_cf + (1-nonmover_pre)*pre_`var'_post_cf
			
			* observed outcome 
			qui egen pre_`var'_obs = sum(`var'*w_pre_obs)
			
		* DP test
	
			* counterfactual
			qui egen dp_`var'_cf = sum(`var'_cf*w_dp_cf)
				
			* observed outcome 
			qui egen dp_`var'_obs = sum(`var'*w_dp_obs)

		}

	format %10.0fc freq* sumfreq*
		
	* build results table
	gen helper = 1
	collapse (mean) nonmover_pre pre_* dp_*, by(helper)

	foreach i in d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat{
		gen diff_`i' = pre_`i'_obs - pre_`i'_cac
		gen diff_pc_`i' = diff_`i' / pre_`i'_cac
		
		gen dp_diff_`i' = dp_`i'_obs - dp_`i'_cf
		gen dp_diff_pc_`i' = dp_diff_`i' / dp_`i'_cf
		}
		
	foreach x in d_age d_male d_ambulance charlindex los2010 n_em_admit2010 d_death_hat d_admit_hat{
	su diff_`x'
	return scalar d_`x' = r(mean)
	su diff_pc_`x'
	return scalar d_pc_`x' = r(mean)
	
	su dp_diff_`x'
	return scalar d_dp_`x' = r(mean)
	su dp_diff_pc_`x'
	return scalar d_pc_dp_`x' = r(mean)
		
		
		}
	
	restore

* END PROGRAMME DEFINE HERE	
end	


#delimit ;
* Simulate 199 times;
simulate diff_male = r(d_d_male) diff_pc_male = r(d_pc_d_male) 				
diff_age = r(d_d_age) diff_pc_age = r(d_pc_d_age)					
diff_ambulance = r(d_d_ambulance) diff_pc_ambulance = r(d_pc_d_ambulance)		
diff_charlindex = r(d_charlindex) diff_pc_charlindex = r(d_pc_charlindex)		
diff_los2010 = r(d_los2010) diff_pc_los2010 = r(d_pc_los2010)	
diff_n_em_admit2010 = r(d_n_em_admit2010) diff_pc_n_em_admit2010 = r(d_pc_n_em_admit2010) 
diff_admit_hat = r(d_d_admit_hat) diff_pc_admit_hat = r(d_pc_d_admit_hat)		
diff_death_hat = r(d_d_death_hat) diff_pc_death_hat = r(d_pc_d_death_hat)		
, reps(199) seed(32786105): myboot;

#delimit ;
* Boostrap;
bstat, stat(observe_d_male, observe_pc_d_male, 		 
observe_d_age, observe_pc_d_age, 			
observe_d_ambulance, observe_pc_d_ambulance,
observe_charlindex, observe_pc_charlindex, 
observe_los2010, observe_pc_los2010, 
observe_n_em_admit2010, observe_pc_n_em_admit2010,	
observe_d_admit_hat, observe_pc_d_admit_hat, 			  
observe_d_death_hat, observe_pc_d_death_hat);


#delimit cr
* Output results
putexcel set "$results/Table2.xlsx", replace
putexcel A1 = "Variable"
putexcel B1 = "Coefficient"
putexcel C1 = "Std. error"

matrix b = e(b)'
putexcel A2 = matrix(b), rownames

matrix se = e(se)'
putexcel C2 = matrix(se)


